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1 Abstract 

In this paper we show by means of numerical experiments that the error intruduced in a 
numerical domain because of a Perfectly Matched Layer or Damping Layer boundary treat- 
ment can be controlled. These experimental demonstrations are for acoustic propagation 
with the Linearized Euler Equations with both uniform and steady jet flows. The propagat- 
ing signal is driven by a time harmonic pressure source. Combinations of Perfectly Matched 
and Damping Layers are used with different damping profiles. These layer and profile com- 
binations allow the relative error introduced by a layer to be kept as small as desired, in 
principle. Tradeoffs between error and cost are explored. 


2 Introduction 

Aeroacoustics research is conducted with theoretical analysis, experimental observations, and 
computation. Computation is being used more often because of improvements in numerical 
methods and computer resources, and is increasingly valuable as a source of insight and 
as a tool for design [3, 18, 23, 28]. Artificial boundary treatments are a critical element 
for Computational AeroAcoustics (CAA), since they are intended to be used to restrict a 
numerical domain whithout causing unacceptable error in the numerical solution. A poor 
boundary treatment can produce errors in the numerical simulation that grow to the size of 
the solution itself, or larger, making the numerical solution of only qualitative or aesthetic 
interest. A poor boundary treatment that might produce small errors could be very costly, 
making an accurate numerical solution practically unobtainable. A good boundary treatment 
is particularly important for simulating time dependant phenomena with detailed physics, 
where experimants may be difficult because of cost or impossible because of scale and detail, 
experiments such as for investigating jet noise. The recent reviews of the extensive work 
done to date on this problem include [2, 5, 6, 13, 15, 16, 24, 27]. In this paper we are 
concerned with Perfectly Matched Layer (PML) and Damping Layer (DL) treatments for 
acoustic propagation by means of Linearized Euler Equations (LEE) in two space dimensions 
with uniform and steady jet flows. A PML boundary treatment is designed for this system by 
ensuring that the plane wave solutions for the LEE in the numerical domain perfectly match 
the plane wave solutions for the adjacent PML domain along the interface between them. 
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Control of the error introduced back into the numerical domain by the PML is attempted by 
a combination of numerical domain filtering and interface conditions, and the PML domain 
damping terms, damping rates and profiles, layer width, and outer boundary treatment. 

In [11] we compared three PML treatments for this problem. Hagstrom [14] has used 
transform methods to analytically formulate a PML boundary treatment (PML A) by oper- 
ator approximation for this system, which has also been further analysed by Motamed [19] 
and Froling [4] . Since PML A has been developed with linear analysis for constant coefficient 
problems, it is not readily apparent how it can be adapted to problems with nonuniform base 
flows, or to fully nonlinear problems. Kami [17] has used eigensystem analysis to develop a 
PML treatment (PML B) that use directional damping in an artificial boundary layer. Kami 
[17] has observed that the damping rate need not be the same for the different eigenfunctions 
of the hyperbolic system, but that the damping terms should not alter the eigenvectors of 
the system with respect to the direction of damping. In [11] we showed that a Thompson like 
characteristic analysis [25, 26] can be used to easily obtain a general damping layer formula- 
tion that does not alter the eigenvectors of the Hyperbolic system. The eigensystem analysis 
that Kami has used could be applied locally, and is not restricted to constant coefficient 
linear problems. The familiar DL or sponge zone treatment does not alter the eigenvectors 
of the system, and the damping profile can be chosen to impose a smooth interface between 
the numerical domain and the sponge zone or damping layer. A DL boundary treatment 
has been recently discussed in [21], where the damping term crU is called a relaxation source 
term, and the damping coefficient a = 1/r, where r is a relaxation time. Bodony [1] has 
also recently given an analysis of sponge zones for computational fluid dynamics. We have 
called this familiar, simple, and easily generalized treatment PML C, even though it is not 
properly a PML. 

In [11] we showed that PML A and PML B were reasonably similar in performance, with 
both producing disturbances that are one to three orders of magnitude smaller than those 
produced by PML C for similar levels of damping. In this paper we consider both PML B 
and C treatments by themselves, and we introduce and test a combination of PML B and 
PML C that is intended to combine the best features of both. This combined PML treatment 
will be called PML BC. PML BC has PML B terms predominate near the interface, and uses 
a smooth transition away from the interface to an outer damping treatment that is purely 
PML C. PML BC is essentially the same as a PML B with different damping profiles for 
each directional eigenvector, and is intended to incorporate the best features of both separate 
PML treatments, the accurate and nondistorting interface of PML B with the nonreflecting 
omnidirectional damping of PML C near the outer PML boundary. Comparisons will be 
made of the relative errors for combinations of damping profiles and layer widths. 


3 PML/DL for the LEE with Uniform Base Flows 

The nondimensionalised Linearized Euler Equations (LEE) in two space dimensions with 
a uniform base flow (ub,Vb) = (M x ,M y ) and a source can be written for the perturbation 
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pressure and velocity U = ( p,u,v) T as 

G, 

‘ My 0 1 \ 

0 My 0 . 

,1 0 My) 

In our numerical experiments we always use the pressure source G = (g, 0,0) T , with 

g(x, y, t ) = 0.01 sin[27rt] exp[— 36(x 2 + y 2 )\. 

This compact radiating source propagates a time harmonic pressure disturbance with max- 
imum values that are O[10 -4 ]. These equations apply in the numerical domain 


dU . dU „ dll 
— + A c — + B c — = 
dt dx dy 

where A c and B c are the constant matrices 


/ 

M x 

1 

° \ 

t 


1 

M x 

° 

, and B c = 

V 

0 

0 

M x ) 

l 


^ n = {(x,y) e [x L ,x R ] x [y B ,y T }}, 

where we are concerned with obtaining accurate results. The numerical domain is surrounded 
with a combination of PML and DL domains, which we take to be Cartesian with linear 
boundaries. The total combined domain is 


& = {<T, y) e \x L - w L , x R + w R \ x [y B - w B , y T + W T ]}, 

where Wl and Wr are the PML/DL domain widths on the left and right, and w B and Wt are 
the widths on the bottom and top. The PML/DL or damping domain is just 

— D — h2jv- 


We will use the homogeneous LEE in the PML/DL domain Q D , since the source term g is 
virtually zero outside of the numerical domain A source in Qjy near the interface with 
Dd could be accomodated merely by including the source in the damping domain equations. 
For this constant coefficient case, a PML B treatment [17, 11] perpendicular to the x axis 


is 


dU A DU „dU 
vyr + — \ - 5 C - — b Scr(x)A c U 

at ox ay 


0, 


where c(x) is the damping profile, 6 the damping scale, and A c the damping coefficient 
matrix. Similarly, a PML B treatment perpendicular to the y axis is 


dU . dU r> dU . . . n 

~m + Ac & + Bc ^ + 


0 . 


Note that since PML B is derived from a directional eigenvector analysis [17, 11], <5 should 
be positive for a layer on the right, and negative for a layer on the left, and similarly for up 
and down. A PML C treatment can be written simply as 


dU dU 
dt ^ c dx 


dU 

B + Sail = 0, 
dy 
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where a would depend on the variable across the layer in either coordinate direction, with 
5 > 0. In any of the fonr corners, the PML/DL treatments from the two adjacent sides can 
be combined by simply adding them. For PML B, this additive corner treatment is just 


OU A 8U n d U . , . . . / \ T> \ TT 

~dt + Ac ~dx + Bc ~dy + S + cr (^) jBc ) U 


0, 


and for PML C it is 


dU A dU n dU . . . . . 

-m +A ^ + B ^ + Ha[x) + a{y))u 


0. 


A more general approach to corners for PML C is to simply make the damping profile 
multidimensional, with 


ou A du „du . . n 

— + A ,— + B c — + Sa(x, V ) u =0. 


The additive corner treatment is then just a simple example of a multidimensional damping 
profile. In [11] we presented an approach for damping in a corner by blending two PML B 
treatments from the sides. In this paper we will use additive corner damping for both PML B 
and C, and the multidimensional corner damping profile for some cases with PML C. There 
appears to be little difference in the errors produced by the various corner treatments, and 
each of these three approaches produces no significant extra error than is produced by the 
PML/DL treatments on the sides. 

A Zero Boundary Data (ZBD) boundary condition is used on all outer boundaries of the 
damping layers, with all data set equal to zero along the entire outer boundary. This is a 
combined Dirichlet and Neumann condition since we keep and propagate spatial derivative 
data at each grid point. Our rationale is that we are primarily concerned with the layer 
treatment, and that if there were a good boundary condition, then the entire layer treatment 
industry would be unnecessary. In addition, if any signal crosses a damping zone and reaches 
an outer boundary with ZBD boundary condition, then the boundary condition will create a 
reflected signal with approximately the same strength. This worst of all backscattering from 
the outer boundary will show how well the layer treatment is performing. This is particularly 
significant for PML B which damps data for outgoing eigenfunctions and amplifies data for 
incoming eigenfunctions. Note that we do not use any filtering to eliminate this effect, but 
just the damping in the formulation of the PML/DL layer treatment. The backscattering 
can be reduced by using a better outer boundary treatment. 

There are various ways to implement a PML or DL. One approach is to simply view 
the damping terms as ordinary terms in the system, and treat them the same as any other. 
For complex algorithms this could lead to significant extra effort. A simple, inexpensive 
and general approach is described by Romenski, Titarev, and Toro in [21]. For the vector 
variable U, consider the evolutionary system 


dU 

~dt 


+ n(c7) = o, 
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where II is the propagator for U, with the associated damped system 

dU 

— + n(u) + su = o, 

where 5 is the damping rate. This is approximated in [21] by 

fjn+l __ fjn ^ 

— = -n (U n )^SU n+1 , or (l + SAt)U n+1 = U n - AtU(U n ). 

Any explicit algorithm for an undamped evolutionary system can be written as 

U n+1 = P(U n ), 

where the algorithm could be to any order in time, with any number of time substeps. The 
damping in [21] can now be adopted as 

(1 + 5At)U n+l = P(U n ), or U n+l = , P ^ U '} , . 
v v ' (1 + 6 At) 

If the damping terms are not simple, say in the form 

dU 

— + U(U) + DU = 0, 

where D is a matrix, such as in PML B, then we can write 

(1 + A tD)U n+l = P(U n ), or U n+l = (1 + A tD)~ l P{U n ). 

The inversion of the matrix damping term can be done locally at each grid point. These 
formulations can be implemented as post processing step after each time step, or after each 
time substep, or set of substeps. The damping terms for PML B and PML C have been 
implemented with algorithmic treatments that are the same as any other terms with all of 
the details required for full accurate time evolution, and also as this simple post processing 
step, and there is little difference in the numerical domain solution, but there can be a large 
difference in the required effort. 

We have generally used four damping profiles. Let £ = x or y, depending on whether or 
not a damping term is for a layer perpendicular to the x or y axis, and let £j be the location 
in x or y of the interface between the numerical domain Vt N and the damping domain f Id- 
We assume that all profiles are set to zero inside the numerical domain. The first damping 
profile was used for the constant coefficient cases in [11], with 

(£-6) 2n 

ffi) (4,w) . 2 . / £ f >2\n’ 

( w 2 + (£ - £/) 2 ) n 

for all £ inside the damping domain, where w is a scaling factor. Note that the order of the 
zero of <7 q at £/ and the smoothness of the interface can be controlled by the choice of n. 
Note also that ctq asymptotes to one, instead of growing indefinitely like typical polynomial 
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damping profiles. The damping profile is multiplied by the maximum damping rate or 
amplitude S. The damping amplitude <5 combines with the layer width scale w to control 
the total damping in the layer. The second damping profile is 

0i(6w) = I w 1 I — I , for |6| < |f| < |6| +W, 

and 

<Ti(Z,w) = !, for 161 + w < |f|- 

Note that cri(£/, w) = 0 , and eri(£j + w,w) = 1 , and that the first 2 n — 1 derivatives of <j\ 
are 0 at both points. The third damping profile is similar to the second, with 

<T2«,«0 = (6^) 2 " ( 2i "~ ( i~ 6> ) 2 '‘- f ° r I 5 '' ^ l«l ^ + 2W . 

and 

cr 2 (£, w) = 0, otherwise. 

Note that cr 2 (6> w ) = 0 = <t 2 (£/ + 2 w,w), and <72(6 + w,w) = 1 , and that the first 2 n — 1 
derivatives of cr 2 are 0 at all three points. Damping profile cr 2 is used for PML B, and <j\ for 
PML C, with n — 3 and additive corner treatments for both profiles. The fourth damping 
profile is 

03(6 w) = exp [-w 2 /£ 2 }, 

where ^ = \x — Xi\ on the right or left, £ = \y — yj\ on the top or bottom, and in the four 
corners £ = yjjx — op) 2 + (y — yj ) 2 . The multidimensional corner treatment is used with (73 
by design. Note that 03 is infinitely smooth at the interface between the numerical domain 
and the damping layers. 

PML BC is a blend of PML B and PML C using a combination of terms with damping 
profiles <7i and <j 2 . For a PML BC perpendicular to the x axis on the right, the damping 
terms are 

D(x)U = S B a 2 (x,w B )A c U + 6 c cri(x, w c )IU, 

where / is the identity matrix. The first term is just PML B for xi < x < xi + 2 w, and 
the second term is PML C. The use of PML B near the interface is intended to ensure that 
lower errors are introduced by the PML, and the final persisting use of PML C is intended to 
provide omnidirectional damping with no amplification of errors back towards the numerical 
domain. A slightly modified form of PML BC delays the introduction of the third term for 
a distance w g from the interface, and can be written as 

D(x)U = S B a 2 (x,w B )A c U + 5 c cri(x - w g , w c )U. 

This is the form that we prefer, with appropriate modifications for right and left, and up 
and down PMLs. Note that there are five parameters, the two damping scales, 5 B and Sc, 
the two width scales w B and wc, and the ”g a p” scaling w g before the third pure PML C 
term becomes active. This PML BC uses a PML B in the layer or zone near the interface 
for accuracy, with greater accuracy from a smaller damping coefficient S B , and PML C away 
from the interface to prevent back scattering from the outer boundary, with greater damping 
from a larger damping coefficient Sc- The interplay of these five parameters will begin to be 
explored in the numerical experiments reported below. 
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4 Numerical Algorithm 

All of the experiments that are reported in this paper are conducted with the Hcrmite / Cauchy- 
Kovalevsky /Taylor method (see [7], [8], [9], [10], [12], [11]). This method uses two staggered 
uniform grids, offset by half a grid spacing in both x and y. The first or basic grid is used 
for initial data and for the numerical solution at each full time step. The second or offset 
grid is used for the numerical solution at each half time step. The time evolution through 
a half time step from one grid to the next starts with multidimensional Hermite spatial 
interpolation of the local data on the current grid stencil, with a tensor interpolant for a 
rectangular grid. The half time step then proceeds by using the governing equations for 
Cauchy-Kovalevsky recursion to produce all of the required time derivatives from the spatial 
derivatives produced by the interpolant. The half time step ends by propagating through a 
half time step with a Taylor series. This approach is not simply an algorithm, but is actually 
a method for developing or specifying numerical algorithms. The recursion routines are de- 
pendant upon the governing equations, but they can be written independantly from the order 
of the method. Various methods can be realized for a particular system by simply swapping 
the interpolation routines, using different stencils and data. This approach has had various 
names, but in this paper, we will call it the Hermite/Cauchy-Kovalevsky/Taylor method, 
or the HCKT method. In this paper we use the c2ol method, which is a HCKT method 
that uses a four point 2x2 grid cell, and interpolates with a bicubic Hcrmite interpolant at 
the cell center, with data for {/, f Xl f yi f xy } at each grid point for any evolving variable /. 
This is a tensor interpolant which simultaneously and consistently calculates all derivatives 
in the local bicubic spatial expansion. The spatial interpolant is used as a local initial data 
surface, and all required time derivatives are computed from the local data surface by means 
of a Cauchy-Kovalevsky recursion. The new solution data is obtained at the cell center 
with a Taylor series in time. The time expansion is to fourth order terms, but because of 
the bicubic tensor interpolant, the overall order is third in space and time, since a fourth 
order time expansion requires fourth order spatial data which the cubic interpolant does not 
provide. Note that the local time evolution is exact up to the order of the method. For two 
dimensional hyperbolic systems, the local data is correctly propagated along characteristic 
surfaces up to the order of the method. In this sense, the HCKT methods can be said to 
correctly extend the method of characteristics to multidimensional hyperbolic systems. By 
using the correct Cauchy-Kovalevsky recursion, this has been done for nonlinear problems 
as well, in particular for convectivcly dominated compressible Navier-Stokes equations. 


5 Results for the LEE with a Uniform Flow 


The numerical experiments in this section are for the LEE with ( pb,Ub,Vb ) = (1,0.4, 0) as 
the base flow, and with the pressure source 


dU . dU „ dU 

dt dx dy 


G, 


where G = ( g , 0, 0) 1 = (0.01 sin[27rt] exp[— 36(x 2 + y 2 )], 0, 0) T . The propagated variables are 
all initially set to zero. All of the results are with the c2ol HCKT algorithm, which is third 
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order accurate in space and time. The algorithm is run with a Taylor series time expansion 
that has fourth order terms, but the third order Hcrmite interpolant restricts the overall 
algorithm to third order. The data that we report is 


IIp-PrIUav = max{| p(x,y,t F ) -p R (x,y,t F )\ : (x,y) G Ojv}, 

where tp is the final time for the simulation being considered, and where p R is a reference 
solution produced with no damping of any kind on a Big Enough Domain (BED) so that 
the solution in is undisturbed at t F by any reflection back from the domain boundary. 


Table 1: \\p — PrIIoo.Ojv a t t = 15, from PML B with oy. 



W B = 1 

w B = 3 

w B = 5 

w B = 7 

10.000 

4.2233 x 10~ 7 

1.5443 x 10~ 7 

3.7913 x 10~ 8 

8.7781 x 10~ 9 

1.000 

8.8739 x 10“ 9 

2.6081 x 10“ 9 

4.5782 x 10“ 10 

1.5450 x 10“ 10 

0.100 

7.2401 x 10” 10 

1.0766 x 10" 10 

2.1359 x lO” 11 

8.5276 x lO" 12 

0.010 

7.3263 x 10” 11 

9.2124 x 10- 12 

2.4369 x 10" 12 

7.8338 x 10” 13 

0.001 

7.3414 x 10“ 12 

9.0569 x 10“ 13 

2.4700 x 10" 13 

7.7644 x 10~ 14 


Note that ||pr||ooAv = °[ 10 4 ]- 


Table 2: \\p — Pi?||oo,fijv a t t = 15, for PML C with oy. 


5c 

w c = 1 

w c = 3 

w c = 5 

w c = 7 

wc = 9 

10.000 

2.5833 x 10“ 5 

1.6333 x 10" 6 

8.2610 x 10~ 7 

6.9061 x 10~ 7 

2.9978 x 10~ 7 

1.000 

2.9492 x 10~ 6 

5.1488 x 10~ 7 

3.5265 x 10~ 7 

1.2282 x 10“ 7 

6.4461 x 10“ 8 

0.100 

2.4761 x 10~ 7 

1.3017 x 10” 7 

4.9640 x 10“ 8 

2.1429 x 10~ 8 

8.7373 x 10~ 9 

0.010 

3.3743 x 10~ 8 

1.3779 x 10~ 8 

5.7895 x nr 9 

2.3194 x 10~ 9 

9.0230 x lO’ 10 

0.001 

3.4936 x HD 9 

1.3810 x nr 9 

5.8907 x 10” 10 

2.3381 x 10~ 10 

9.0522 x 10" 11 


Note that ||pr||ooAv = O[10 4 ]. 


The first series of computations are with PML B and PML C by themselves, and are 
presented in Tables 1 and 2. These calculations are for the numerical domain 

Da j = [—3,7] x [—5,5], with t F = 15, 

with algorithm parameters Ax = 1/24 = Ay, and At = 1/48. The damping domain layer 
widths are w F = w R = w R = w F = 10, the oy damping profile and additive corner treatments 
are used for both PML B and PML C, and the damping terms are fully propagated. Note 
that the propagation speeds of the wave fronts are 1.4 in the +x direction, 0.6 in the —x 
direction, and 1 in the Ay directions, and consequently, that the propagating signal never 
reaches the outer boundary of Vt D by t F = 15. The data in Tables 1 and 2 clearly shows that 
errors introduced by a PML/DL can be controlled, even with a PML/DL on all sides and 
corners of the numerical domain. The three controls that effect the results in Tables 1 and 2 
are the type of layer that is used, the damping amplitude d, and the damping width scale w. 
For every combination of d and w , PML B produces errors that are between two and three 
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orders of magnitude smaller than those produced by PML C. The effect of both the damping 
amplitude and width is linear for both PML B and PML C. The solution scale is O[10 -4 ], 
so that the relative error is from O[10 _1 ], or ten percent, to O[10~ 10 ]. Any error that is 
relatively less than 0[1CT 3 ] is invisible at the scale of the solution, and this can be readily 
achieved. It should be noted that if the damping layer widths wb or wc are greater than 5, 
then the effect of the maximum damping profile amplitude is not returned to the numerical 
domain Q/v where ||p — PrHoo.Qjv i s computed. The damping profile can cause an error by 
its initial slope, with the layer scale and profile type both playing a role, by the maximum 
damping amplitude, and possibly by higher derivatives from the profile shape. The smallest 
errors in Tables 1 and 2 are produced by combinations of small damping amplitude and 
large damping width scales. The total damping is a product of the amplitude and width, 
so that a small damping amplitude needs a corespondingly wide damping layer to provide a 
given total damping. The very smallest errors are produced by the very smallest damping 
amplitudes, and they provide very little damping except over very wide damping layers. 


Table 3: \\p-pR\\oo,n 

N at t — 25, for PML BC, with 

(j 2 and wc = 3. 

w c 

wr — 1 

w B = 3 

w B = 5 

w B = 7 

1 

1.2269 x 10~ 6 

1.2225 x 10“ 6 

1.2278 x 10~ 6 

1.2290 x 10" 6 

3 

4.1117 x 10“ 7 

4.1009 x 10" 7 

4.0892 x 10“ 7 

4.1163 x 10~ 7 

5 

1.4603 x 10~ 7 

1.4482 x 10" 7 

1.4424 x 10“ 7 

1.4563 x 10“ 7 

7 

4.4231 x 10~ 8 

4.2385 x 10" 8 

4.2241 x 10~ 8 

4.3100 x 10“ 8 

9 

2.0671 x 10“ 8 

1.9703 x 10" 8 

1.9697 x 10“ 8 

1.9843 x 10~ 8 


Note that ||pi?||oo,n 

3 

o 

o 

i 


Table 4: ||p - P/?||oo,n 

^ at t — 25, for PML BC, with 

ex i and wc = 5. 

w c 

w B = 1 

w B = 3 

w B = 5 

w B = 7 

1 

6.8972 x 10“ 7 

6.8669 x 10" 7 

6.8458 x 10“ 7 

6.8825 x 10~ 7 

3 

7.2773 x 10~ 8 

7.1685 x 10" 8 

7.1769 x 10~ 8 

7.2137 x 10~ 8 

5 

1.6277 x 10“ 8 

1.6506 x 10“ 8 

1.6558 x 10“ 8 

1.6707 x 10~ 8 

7 

5.4811 x nr 9 

4.0522 x HC 9 

4.0495 x HC 9 

4.2073 x HC 9 

9 

5.0071 x 10“ 9 

1.5652 x KC 9 

1.2252 x nr 9 

1.3843 x HC 9 


Note that ||pr||ooAv = °[ 10 4 ]- 


The second series of computations uses the combined PML BC damping layer, with 
results presented in Tables 3 and 4. These calculations are for the numerical domain 

Q/v = [-6, 14] x [-10, 10], with t F = 25. 

The algorithm parameters are Ax = 1/24 = Ay, and At = 1/48, the damping domain 
widths are wl = wr = wr = w F = 10, the damping terms are fully propagated, and an 
additive corner treatment is used. The propagating signal does not return from the outer 
boundary of the damping domain by the final simulation time t F = 25. Recall that a- 2 is a 
pulse, so that PML B is only active close to the boundary, and that PML C applies only at 
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a distance greater than wq from the damping domain interface. The damping amplitudes 
5b = 1 for the PML B terms and 5c — 1 for the PML C terms have been chosen as moderate 
but relatively effective damping amplitudes. Note that the data in both tables shows that 
the PML BC combination is insensitive to the length scale of the PML B terms. Note also 
that the PML C terms decrease the effectivness of the PML B terms, and that the PML BC 
errors are from one to three orders of magnitude larger than the comparable errors for PML 
B by itself. Comparing the results in Table 3 with Table 2 shows that these PML BC results 
are slightly better than the PML C results by a factor of about three, and that the PML 
BC results are worse than the PML B results by from one to three orders of magnitude. 
This suggests that the PML BC error in these computations is dominated by the effect of 
the PML C terms. Comparing the results in Table 4 with Table 1 shows that these PML 
BC results are worse than the PML B results by about an order of magnitude, and that the 
PML BC results are better than the PML C results by about an order of magnitude. The 
best performance of PML BC has PML C fully engaged at a distance of w g +wc = 12 or 14 
from the interface. If PML C had been used by itself with w g = 0 and wc = 12 or 14, then 
the relative advantage of PML BC would have been significantly reduced. The PML BC 
combination can provide some damping layer performance improvement over a pure PML 
C, but it appears from these initial results that this combination will work well only with a 
relatively large gap before the PML C terms take effect. The best strategy for employing the 
PML BC combination seems to be to reduce the solution amplitude by about three orders 
of magnitude with the PML B terms, and then to use the PML C terms only after most or 
all of this initial damping has been accomplished. The PML BC combination will require 
wide damping layers to be effective. 


Table 5: \\p — PrIIoo.Ujv at t — 35, for PML C with PML width 15 and cr 3 . 


<5c 

W C = 1 

w c = 5 

wc = 9 

wc = 13 

10.00 

8.6856 x 10” 6 

4.4231 x 10" 7 

2.6493 x 10~ 7 

1.6355 x 10~ 7 

1.000 

5.4827 x 10~ 7 

1.4929 x 10“ 7 

7.8054 x 10“ 8 

3.2457 x 10“ 8 

0.100 

6.9008 x 10~ 8 

3.0104 x 10~ 8 

1.3658 x 10~ 8 

6.6209 x HR 9 

0.010 

1.0729 x 10~ 8 

4.0779 x 10~ 9 

1.7786 x 10~ 9 

7.3515 x 10“ 10 

0.001 

1.1433 x 10~ 9 

4.2478 x 10~ 10 

1.8288 x 10~ 10 

7.4310 x 10' 11 


Note that ||pr||ooAv = °[ 10 4 ]- 


The third series of computations is with PML C by itself, and is presented in Table 5. 
These computations use the a 3 damping profile, the multidimensional corner treatment, and 
the Romenski, Titarev, and Toro [21] implementation of the damping terms. The numerical 
domain is 

CIn = [—6, 14] x [—10, 10], with tF = 35, 

and the damping domain widths are wl = wr = wb = u>t = 15. The algorithm parameters 
are Ax = 1/12 = Ay, and At = 1/24. The Damping zone is large enough so that the 
propagated signal does not return by t F = 35 from the outer boundary of the damping 
domain, so that the outer boundary does not effect the solution in the numerical domain. 
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The data in Table 5 is comparable to the results in Tables 1 and 2, in the sense that 
all three tables show the effect of the damping treatment without any reflection from the 
outer boundary of the damping layer. Note that the damping profile scales in Table 2 are 
wc = 1, 2, 3, 5, and 7, but in Table 5 they are wc = 1, 5, 9, and 13. It is immediately clear that 
the results in Table 5 are only slightly different from the comparable results with the same 
damping profile scale for PML C in Table 2, with errors that are smaller by no more than 
a factor of two. The slightly better results in Table 5 could possibly be due to the greater 
decay of the wave front from the initial impulsive start because of the longer simulation 
time. On the other hand, for any particular damping amplitude Sc, damping profile a 3 
increases about three times more gradually than a i, so that the a 3 profile scale is effectively 
three times larger than Wq- This more gradual increase in the damping profile will tend to 
produce somewhat smaller errors in the numerical domain. Nonetheless, the similarity of 
the results in Tables 2 and 5 shows that the damping profile, grid resolution, domain size 
and simulation time, corner treatment, and implementation of the damping terms all have 
no significant effect upon the accuracy of the damping treatment. Note in this regard, that 
all of the damping profiles have been chosen so that, at the interface between the numerical 
and damping domains, all of the spatial derivatives of the profiles are zero up to at least the 
order of the numerical method. In particular, a 3 is C°° with every spatial derivative equal to 
zero at the interface, so that the interface between the numerical domain and the damping 
layer provides a smooth transition between the governing systems in each domain. One of 
the implications of the data that has been presented so far is that relatively small errors 
are introduced into the numerical domain solution if the transition from the numerical to 
the damping zone is gradual as well as smooth. The data in Table 5 also serves to validate 
the algorithm mix that is used here, and to show that the algorithm mix and simulation 
parameter selection is representative of the other results that have been discussed. 

The fourth series of computations is with PML C by itself, and is presented in Tables 6, 
7, 8 and 9. These computations use the same algorithm mix that produced Table 5, with 
the same damping profile, multidimensional corner treatment, and implementation of the 
damping terms. The numerical domain is the same, 

n N = [-6,14] x [-10,10], 

with the same algorithm parameters, Ax = 1/12 = Ay and At = 1/24. However, the final 
simulation times and the damping layer widths are all different. Table 6 is for final simulation 
time t,p = 55 with damping domain widths wl = wr = wb = wp = 15, Table 7 is for tp = 65 
with wl = wr = wb = u)p = 20, Table 8 is for tp = 75 with wl = wr = wb = wp = 25, and 
Table 9 is for tp = 85 with wl = wr = wb = wp = 30. In each case, the final simulation 
time is large enough so that the transient effects from the impulsive start can pass out of the 
numerical domain, and so that the propagated signal can reflect from the outer boundary 
of the damping domain back through the entire numerical domain. The results in Tables 
6-9 are designed to show the effect of reflection from the outer boundary as mediated by the 
damping profile and the damping zone width. For almost all combination of Sq and Wc in 
Tables 6-9, \\p— PkIIoo.Hjv decreases as the damping layer width increases, with exceptions that 
appear to be minor variations around a limiting value. In these four tables, insufficient total 
damping for a fixed Sc is indicated if ||p — p/j||oo,n w increases with wc, which stretches the 
darnpng profile and therefore decreases the total damping in the layer. Similarly, insufficient 
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Table 6 : \\p — PR||oo,fijv at t — 55, for PML C with PML width 15 and < 73 . 


5c 

w c = 1 

w c = 5 

wc = 9 

wc = 13 

10.0 

8.4234 x 10 " 6 

9.9517 x 10 “ 8 

5.9050 x 10 ~ 8 

9.2317 x 10 ~ 8 

1.0 

5.3772 x 10 ” 7 

3.8008 x 10 ~ 8 

5.0535 x 10 “ 8 

6.8419 x 10 ~ 8 

0.1 

1.9056 x 10 “ 7 

6.2012 x 10" 7 

1.3412 x 10 ' 6 

2.1143 x 10 ' 6 


Note that ||pr||ooAv = °[ 10 4 ]- 


Table 7: \\p — PrIIoo,^ at t — 65, for PML C with PML width 20 and 0 - 3 . 



w c = 1 

w c = 5 

wc = 9 

wc = 13 

10.0 

8.3917 x 10 “ 6 

1.6884 x 10 ~ 8 

2.2645 x 10 “ 8 

3.7814 x 10 “ 8 

1.0 

5.4330 x 10 ~ 7 

1.5659 x 10 “ 8 

2.6035 x 10 “ 8 

2.3734 x 10 ~ 8 

0.1 

7.2112 x 10 ~ 8 

1.6864 x 10 ~ 7 

4.3244 x 10 “ 7 

8.3129 x 10 " 7 


Note that ||pr||ooAv = ^[lO 4 ]- 


Table 8 : \\p — p^Hoo,^ at t = 75, for PML C with PML width 25 and 03 . 


5c 

w c = 1 

w c = 5 

wc = 9 

wc = 13 

10.0 

8.3890 x 10 “ 6 

2.0426 x 10 “ 8 

1.0227 x 10 ~ 8 

8.0675 x 10 " 9 

1.0 

5.4134 x 10 ~ 7 

9.5706 x 10 " 9 

6.1355 x 10 “ 9 

1.1910 x 10 ~ 8 

0.1 

4.7353 x 10 ~ 8 

4.7189 x 10 ~ 8 

1.3228 x 10 “ 7 

2.9477 x 10 ~ 7 


Note that ||pr||ooAv = °[ 10 4 ]- 


Table 9: \\p — p^Hoo,^ at t — 85, for PML C with PML width 30 and cr 3 . 


5c 

w c = 1 

w c = 5 

wc = 9 

wc = 13 

10.0 

8.3996 x 10 ~ 6 

2.9749 x 10 “ 8 

1.2372 x 10 “ 8 

8.4349 x 10 “ 9 

1.0 

5.4403 x 10 ” 7 

1.3398 x 10 “ 8 

6.3281 x 10 “ 9 

3.4250 x 10 " 9 

0.1 

4.6392 x 10 ~ 8 

1.5651 x 10 “ 8 

4.0343 x 10 “ 8 

9.9061 x 10 “ 8 


Note that ||pr||ooAv = °[ 10 4 ]- 


total damping for a hxed wc is indicated if \\p — p/jHoo.njv increases with a decrease in Sc, 
which diminishes the damping profile, decreasing the total damping. Note from the data for 
wq = 1 in Tables 6-9 that the errors for Sc = 1 and 10 are the same for any PML width, 
from 15 to 30, and that these errors are close to the comparable errors in Table 5 where 
the outer boundary has no effect. We may conclude that a PML width of 15 is sufficient 
to produce the best result possible with wc = 1 and Sc = 1 or 10 , or that the profile itself 
creates an error level that is not diminished by further damping farther from the interface. 
By the time the solution has propagated through a layer of width at least fifteen, it has 
been damped sufficiently so that whatever error is reflected from the outer boundary is less 
than the error produced by the damping layer with these profiles. Similarly, for wc = 1 and 
Sc = 0 . 1 , a damping layer width of about 20 or slightly more is sufficient to produce the best 
that is possible with this profile. Table 10 presents crude estimates of damping layer widths 
that suffice to produce the best possible results with the damping profile and parameter 
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combinations that were used to produce the error data in Tables 8-9. A total damping can 


Table 10: Layer width for smallest errors with PML C and <73. 


<5c 

w c = 1 

w c = 5 

wc = 9 

wc = 13 

10.00 

15 

20 

25 

25 

1.000 

15 

25 

25 

> 30 

0.100 

20 

30 

> 30 

> 30 


be defined as 

TD(w d ,S c ,w c ) = S c 

If the results in Tables 6-9 are sorted into smaller and greater errors than the comparable 
cases in Table 5, then 

TD(wc,wd ) > 1.7 for smaller errors, 


r w D 


<?3{£,wc)d£. 


and 


TD(wc,Wd ) < 1.4 for greater errors. 


This in principle permits at least the crude estimation of damping layer parameters that 
would maintain a specified error bound. It must be remembered that relatively small errors 
are produced if the transition from the numerical domain to the damping zone is gradual as 
well as smooth, or with large damping profile scaling wq and small damping amplitudes Sc- 


6 DL for the LEE with a Nonuniform Base Flow 


The LEE in two space dimensions with a nonuniform base flow requires two thermodynamic 
variables, and we choose to use the disturbance pressure p and specific volume v s = 1/p, 
where p is the density. In terms of the perturbation quantities 


U = (p,u,v,v s ) T , 


the LEE can be written as: the pressure equation 


dp _dp _d p _ du dv dp dp du dv 

ln +U Tx +l Ty + 1P{ lH + Try ) + U l^ + l % + rP{ irx + Try ) 


0; 


the velocity equations 


du _du ~du 1 _ dp du du 1 dp 

H + U d~x +V 'Fy + ywf + + V ~Fy + ysf 


0 . 


and 


dv _ dv _dv 1 dp dv dv 1 dp 

— b U— b V— 1 b U— b V— 1 7rV s — 

dt dx dy 7 Mj. dy dx dy 7MJ. dy 


0; 
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and the equation for the specific volume 


d v s _ dv s _ dv s 

af + u dF + , % 


, du dv 


dv 


S 


dv 


M + di ) + U ~d^ +V dy 


s 


,du 


M, + 


dv ^ 

d ? 


= 0 ; 


where the base flow is (p, u, v, v s ) T , and where Mr is a reference Mach number. The dimen- 
sionless variables for the base flow should satisfy 


a 2 M 2 R =p/p = pv s , 


so that if we want the nondimensional speed of sound to be a = 1, then we must have 


For this paper we assume a constant base temperature, and choose 

p = 1, and v s = 1, so that Mr = 1. 

In this case, we shall write the LEE with nonuniform base flows as 


dU 

~dt 


where A v , B v , and C v are 


A„ — 


( ® 

v sh 

0 

V o 


IV 

u 

0 

-v* 


. dll 

_ dll 


A„—~ + B„— + CM = 0. 

dx 

dy 


0 0 \ 


( v 0 yp 0 \ 

0 0 


0 v 0 0 

u 0 

, B v = 

hs/y 0 v 0 

0 u ) 


V 0 0 -v s v ) 


and 


/ ry ( du I dv\ dp dp 

I ’ \ dx ' dy ) dx dy 

0 


c v = 


V 


0 

0 


du du 
dx dy 

dv dv 
dx dy 


1 dp 
7 dx 

1 dp 
7 dy 


dv s dv s ( du I dv\ I 

dx dy \ dx ' dy ) / 


and where U = (p,u,v,v s ) T are the unsteady perturbation quantities. Note that these 
coefficient matrices could be nonconstant in time as well as in space, and that a minor 
modification would permit a nonuniform base temperature. We considered the LEE with 
the pressure source 

g(x,y,t ) = 0.01sin[27rt] exp[— 36(x 2 + y 2 )}, 

in the pressure equation. With this pressure source, the LEE with nonuniform flow can be 
written as 

9 V A,— + bSS + C V U = G, 


dt 


dx 


dy 
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where G = (g, 0, 0, 0) T . 

In the case of a nonuniform base flow, a PML B treatment perpendicular to the x axis 
will be adapted in the form 

dU OU OU 

— — + Ay— h By— b CyU + 5 (T (x) AyU = 0, 

at ox Oy 

where cr(x) is a variable damping prohle, 5 is a damping scale, and the damping coefficient 
matrix is just the x propagation matrix A v . A similar formulation can be written down for 
a PML B treatment perpendicular to the y axis, 


— — + Ay— b By— b CyU + S <7 (?/) B yU ~ 0. 

at ox oy 

Note that since PML B is derived from a directional eigenvector analysis, 6 should be positive 
for a layer on the right, and negative for a layer on the left, and similarly for up and down. 
A PML C treatment can be written simply as 


dU 

~dt 


dU dU 

Ay— b By— b CyU + 5 <7 U = 0, 

Ox Oy 


where a would depend on the variable across the layer in either coordinate direction, and 

5 > 0. 


7 Results for the LEE with a Parallel Jet 

The next set of experiments are for the LEE with a nonuniform base flow and a pressure 
source. The pressure source is the same time harmonic oscillation of a narrow Gaussian 
spatial prohle that was used above, 

g(x, y, t ) = 0.01 sin[27rt] exp[— 36(x 2 + y 2 )\. 

The nonuniform base how that we use here is a simple parallel jet how, with 

U = (p, ■u,h,h s ) 2 = (1, 0.4 + 0.4exp[— 36y 2 ], 0, 1). 


The numerical domain is 


Ov = {(z,3/)e [-3,7] x [-5,5]}, 

and the third order c2ol HCKT algorithm is used, with Ax = 1/24 = Ay, and At = 1/60. 
Solutions are computed for t < 10. A PML C is used with damping prohle cr 1 and zone width 
Wr = wl = Wb = Wt = 5. The damping zone uses the PML C formulation for nonuniform 
hows in the inflow and outflow damping zones 

C D r = {(x, y) e [-8, -3] x [-5, 5]} and 0 DO = {(x, y) e [7, 12] x [-5, 5]}, 
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Table 10: ||p — Pi?||oo,n w at t — 10, for PML C with o\ and PML width 5. 



W C = 1 

w c = 3 

w c = 5 

w c = 7 

10.000 

1.7287 x 10“ 5 

2.4439 x 10" 6 

4.1319 x 10~ 7 

1.2125 x 10" 7 

1.000 

2.4723 x 10~ 6 

4.1704 x 10” 7 

6.4439 x 10“ 8 

1.5573 x 10“ 8 

0.100 

3.2806 x 10” 7 

4.5465 x 10“ 8 

7.2817 x 10" 9 

1.5973 x 10" 9 

0.010 

3.5859 x 10~ 8 

4.5992 x 10” 9 

7.3718 x lO" 10 

7.3718 x 10" 10 

0.001 

3.6219 x 10“ 9 

4.6049 x lO" 10 

7.3809 x IQ’ 11 

1.6018 x 10" 11 


Note that ||pr||ooAv = O[10 4 ], 


and the PML C formulation for uniform flows in the top and bottom damping zones 

fW = {(x, y ) e [-8, 12] x [5, 10]} and Q DB = {(x, y) e [-8, 12] x [-10, -5]}. 

The corners are treated additively and as if both side zones used the uniform flow PML C 
formulation. The data that will be reported is the maximum absolute error in the numerical 
domain ||p — PR||oo,n w , where p R is a reference numerical solution at t = 10. The reference 
pressure solution p R is O[10 -4 ] in Q?v at tp = 10. 

The data in Table 10 for PML C with a nonuniform flow is comparable to the data in 
Table 2 for PML C with a uniform flow. Both sets of data are for the same numerical domain, 
and the same spatial grid resolution. Because of the higher maximum velocity of the parallel 
jet the time resolution has been increased from At = 1/48 to At = 1/60. Note that this jet 
is superimposed upon the same uniform flow as was used for the experiments that produced 
Table 2. Recall that the data in Table 2 is with a damping zone width = 10 and final 
time tp = 15, while the data in Table 10 is with a damping zone width wp = 5 and final 
time tp — 10. The same damping profile is used in both sets of calculations, with the same 
space scales for the profiles, and the same set of damping amplitudes (except that wc = 9 
has not been used for the parallel jet case). The damping terms have been implemented 
with full time evolution, and not with the simpler postprocessing algorithm, just as for the 
simulations that produced the data in Table 2. We note that almost all of the data for 
Hp-PrIIooAv i n Table 10 is smaller than the corresponding data in Table 2 by up to almost 
an order of magnitude. As in Table 2, the errors in Table 10 range in absolute terms from 
O[10~ 5 ] to O[10 -11 ], or in relative terms from O[10 _1 ] to O[10” 7 ]. Note here as well that the 
PML error decreases linearly with the damping amplitude, and linearly with an increase in 
the damping profile width scale. As in the uniform flow case, the error induced by a PML C 
damping zone with a nonuniform flow is controllable. Overall, as demonstrated by the data 
in Table 10, the performance of PML C with a nonuniform flow is remarkably close to the 
performance of PML C with a uniform flow. 


8 Conclusions 

Computations for the Linearized Euler Equations (LEE) have been reported for acoustic 
propagation with a uniform flow and a parallel jet. Two damping layer treatments for the 
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Linearized Euler Equations have been considered, Perfectly Matched Layer (PML) PML B 
based upon directional damping [17, 11], and Damping Layer (DL) PML C [21, 11], The 
multilayer PML BC has been introduced, with a damping matrix that combines PML B 
for accuracy near a PML interface with PML C for damping in an outer layer. A series of 
numerical experiments have been conducted with all three boundary layers, primarily for 
the LEE with a uniform base flow, but also with a jet flow. The intent of these experiments 
has been to try and understand what are the critical parameters that effectively control 
the errors produced by a PML/DL, and what is the practically achievable level of error 
that can be obtained. Various combinations of numerical domain size, simulation time, 
grid resolution, damping layer size, treatments for layer corners, implementation of damping 
terms, damping profiles types, and damping profile parameter mix have been used in the 
numerical experiments. The data from the numerical experiments shows that: 

1. The errors introduced into a numerical solution by either a PML B or a PML C can 
be controlled. Results from the numerical experiments show relative error reductions 
by as little as OflCD 1 ] to as much as O[10~ 10 ]. It seems clear that smaller errors can 
be obtained with sufficient effort. 

2. Control of the error from a PML/DL appears to be effected most strongly by the layer 
width, and the amplitude and width scale of the damping profile. If the transition from 
the numerical domain to the damping layer is sufficiently smooth, then effect of both 
the damping amplitude and width is linear. Differences in the damping prohle, grid 
resolution, domain size and simulation time, corner treatment, and implementation 
of the damping terms all have much less effect upon the accuracy of the damping 
treatment. If the damping treatment offers insufficient total damping, then reflection 
from the outer boundary of the damping layer can become the dominant source of the 
error produced by the layer. 

3. PML B produces approximately two orders of magnitude less disturbance in the numer- 
ical domain than PML C, but PML B amplifies distortions back towards the numerical 
domain, while PML C dampens omnidirectionally. 

4. The best strategy for employing the PML BC combination seems to be to reduce the 
solution amplitude by about three orders of magnitude with the PML B terms, and 
then to use the PML C terms only after most or all of this initial damping has been 
accomplished. Because of this, the PML BC combination appears to require wide 
damping layers to be effective. 

5. Estimation of the total damping as a function of layer width, damping amplitude, 
and damping spatial scale permit estimation of layer parameters needed to maintain 
a specified error bound. Total damping is a product of the amplitude and width, so 
that a small damping amplitude requires a corespondingly wide damping layer. 

Relatively small errors are introduced into the numerical domain solution if the transition 
from the numerical to the damping zone is gradual as well as smooth. Accurate damping 
layers with O[10 5 ] relative errors are not difficult to achieve. Very accurate damping layers are 
difficult to achieve and appear to require wide layers and substantial computational effort. 
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The treatment of the outer boundary of the damping layer and its effect on controlling 
the error produced by the layer as a whole has not been considered here. Every order of 
magnitude decrease in the error produced by the outer boundary treatment will have a 
significant effect on the error produced by the damping layer, or on the cost of producing 
any specified error level. 
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